本程序读取GRAPES_1km数值模式预报数据, 叠加在三维地形上进行可视化显示.
载入必备的程序库用于读取, 处理和可视化数据.
library(ggplot2)
Use suppressPackageStartupMessages() to eliminate package startup messages
library(raster)
library(png)
library(rayshader)
library(leaflet)
library(metR)
Registered S3 method overwritten by 'data.table':
method from
print.data.table
library(geosphere)
library(magrittr)
载入程辑包:‘magrittr’
The following object is masked from ‘package:raster’:
extract
library(nmcMetIO)
work_dir <- "/media/kan-dai/workspace/workcode/personal/dk_winter_olympic_3D"
setwd(work_dir)
地形数据从"https://www.gmrt.org/services/gridserverinfo.php#!/services/getGMRTGrid" 下载, 输出为geotiff格式, 分辨率选择最大.
# 载入地形数据
elev_file <- "../data/stadium_topo_small.tif"
elev_img <- raster::raster(elev_file)
# 设置坐标信息
crs(elev_img) <- CRS('+init=EPSG:4326')
# 利用leaflet可视化地形数据
pal <- leaflet::colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), # span of values in palette
values(elev_img), #range of values
na.color = "transparent")
leaflet() %>%
addProviderTiles("OpenTopoMap", group = "OSM") %>%
addRasterImage(elev_img,
color=pal,
opacity = 0.6) %>%
leaflet::addLegend(values = values(elev_img),
pal = pal,
title = "Digital elevation model")
# convert it to a matrix which rayshader can handle.
elev_matrix <- raster_to_matrix(elev_img)
[1] "Dimensions of matrix are: 1058x906."
# fill NA values
elev_matrix[!is.finite(elev_matrix)] <- min(elev_matrix, na.rm=TRUE)
计算地形阴影层, 存入变量, 可以在三维地形显示时叠加到地形上. 将地形数据和阴影层叠加起来, 展示二维效果.
shademat_file <- file.path(work_dir, 'data', "shademat_small.rds")
ambmat_file <- file.path(work_dir, 'data', "ambmat_small.rds")
raymat_file <- file.path(work_dir, 'data', "raymat_small.rds")
if(!file.exists(shademat_file)){
shademat <- sphere_shade(elev_matrix, texture = "imhof4")
saveRDS(shademat, file=shademat_file)}
shademat <- readRDS(shademat_file)
if(!file.exists(ambmat_file)){
ambmat <- ambient_shade(elev_matrix, zscale=30)
saveRDS(ambmat, file=ambmat_file)}
ambmat <- readRDS(ambmat_file)
if(!file.exists(raymat_file)){
raymat <- ray_shade(elev_matrix, zscale=30, lambert = TRUE)
saveRDS(raymat, file=raymat_file)}
raymat <- readRDS(raymat_file)
# plot 2D
shademat %>%
add_shadow(raymat, max_darken = 0.5) %>%
add_shadow(ambmat, max_darken = 0.5) %>%
plot_map()
# 获得地图范围
bbox = extent(elev_img)
bbox = list(p1 = list(long = bbox[1], lat = bbox[3]),
p2 = list(long = bbox[2], lat = bbox[4]))
# get map image size
image_size <- define_image_size(bbox, major_dim = 2000)
overlay_file <- "../data/stadium_map_small.png"
if(!file.exists(overlay_file)){
get_arcgis_map_image(bbox, map_type = "World_Topo_Map", file = overlay_file,
width = image_size$width, height = image_size$height,
sr_bbox=4326)
}
载入需要的程辑包:httr
载入需要的程辑包:glue
载入程辑包:‘glue’
The following object is masked from ‘package:raster’:
trim
载入需要的程辑包:jsonlite
{
"results": [
{
"paramName": "Output_File",
"dataType": "GPDataFile",
"value": {
"url": "https://utility.arcgisonline.com/arcgis/rest/directories/arcgisoutput/Utilities/PrintingTools_GPServer/x_____xsJm_6ylOoAUhl7bzx4mC4w..x_____x_ags_cb2a025a-7e61-11ec-b522-22000a88f94a.png"
}
}
],
"messages": []
}
image saved to file: ../data/stadium_map_small.png
overlay_img <- png::readPNG(overlay_file)
plot.new()
# 显示png文件
rasterImage(overlay_img, 0,0,1,1)
# load forecast
#datafile <- "/media/winter_olympic_data/cma_meso_1km/202201220800/t2m/202201220800.001"
#data <- read_micaps_4(datafile, outList=FALSE)
dataT <- retrieve_micaps_model_grid("GRAPES_1KM/TMP/2M_ABOVE_GROUND/", filter="*.001")
# extract region
dataT <- dataT[lon >= bbox$p1$long & lon <= bbox$p2$long & lat >= bbox$p1$lat & lat <= bbox$p2$lat]
# smooth the data
out <- smooth2d(dataT$lon, dataT$lat, dataT$var1, theta=0.02, nx=128, ny=128)
# plot contour image
colors <- c("#3D0239","#FA00FC","#090079","#5E9DF8","#2E5E7F",
"#06F9FB","#0BF40B","#006103","#FAFB07","#D50404","#5A0303")
p_temp <- ggplot(out, aes(x, y, z=z)) +
# geom_contour_fill(breaks=seq(-12, 10, by=0.2)) +
geom_contour(breaks=seq(-12, 10, by=0.2), color="black", size=0.5) +
geom_text_contour(stroke = 0.3, size=16, min.size=3) +
scale_fill_gradientn(name="Temp", colours=colors, limits=c(-12, 5)) +
scale_x_continuous(expand=c(0,0)) +
scale_y_continuous(expand=c(0,0)) +
guides(fill=FALSE)+
theme(panel.spacing=grid::unit(0, "mm"),
plot.margin=grid::unit(rep(-1.25,4),"lines"),
panel.background = element_rect(fill = "transparent"), # bg of the panel
plot.background = element_rect(fill = "transparent", color = NA), # bg of the plot
axis.text=element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
legend.key.size = unit(5, "cm"),
legend.key.height = unit(2.5, 'cm'),
legend.key.width = unit(6, 'cm'),
legend.position = c(0.25, 0.06),
legend.direction="horizontal",
legend.text = element_text(size=32),
legend.title = element_text(size=32),
legend.background = element_rect(fill="lightblue",
size=0.5, linetype="solid"))
# export png file
overlay_file_temp <- "../data/fine_gridded_forecast_temp.png"
png(overlay_file_temp, width=image_size$width, height=image_size$height, units="px", bg="transparent")
p_temp
dev.off()
null device
1
img <- png::readPNG(overlay_file_temp)
plot.new()
rasterImage(img, 0,0,1,1)
colors <- c("#3D0239","#FA00FC","#090079","#5E9DF8","#2E5E7F",
"#06F9FB","#0BF40B","#006103","#FAFB07","#D50404","#5A0303")
p_temp <- ggplot(out, aes(x, y, z=z)) +
geom_contour_fill(breaks=seq(-16, 8, by=0.25)) +
geom_contour(breaks=seq(-16, 8, by=0.25), color="black", size=0.5) +
scale_fill_gradientn(name="Temp", colours=colors, limits=c(-12, 5)) +
guides() +
theme(legend.key.size = unit(1, "cm"),
legend.key.height = unit(0.4, 'cm'),
legend.key.width = unit(1, 'cm'),
legend.position = c(0.1, 0.1),
legend.direction="horizontal",
legend.background = element_rect(fill="lightblue",
size=0.5, linetype="solid"))
p_temp
# load fine gridded forecast
dataU <- retrieve_micaps_model_grid("GRAPES_1KM/UGRD/10M_ABOVE_GROUND/", filter="*.001")
dataV <- retrieve_micaps_model_grid("GRAPES_1KM/VGRD/10M_ABOVE_GROUND/", filter="*.001")
# extract region
dataU <- dataU[lon >= bbox$p1$long & lon <= bbox$p2$long & lat >= bbox$p1$lat & lat <= bbox$p2$lat]
dataV <- dataV[lon >= bbox$p1$long & lon <= bbox$p2$long & lat >= bbox$p1$lat & lat <= bbox$p2$lat]
dataUV <- merge(dataU, dataV, by=c("lon", "lat", "lev", "time", "initTime", "fhour"))
# plot wind streamlines
p_wind <- ggplot(dataUV, aes(lon, lat)) +
geom_streamline(aes(dx=dlon(dataUV$var1.x,lat), dy=dlat(dataUV$var1.y), size=..step..,
alpha=..step.., color=sqrt(..dx..^2 + ..dy..^2)),
L=0.03, res=2, n=10, S=5, lineend="round", size=5, arrow=NULL) +
geom_arrow(aes(dx = dataUV$var1.x*0.8, dy = dataUV$var1.y*0.8),
skip = 1, color = "black", arrow.length=3, size = 2.5) +
viridis::scale_color_viridis(guide="none") +
scale_size(range=c(0,3), guide="none") +
scale_alpha(guide="none") +
scale_x_continuous(expand=c(0,0)) +
scale_y_continuous(expand=c(0,0)) +
theme(panel.spacing=grid::unit(0, "mm"),
plot.margin=grid::unit(rep(-1.25,4),"lines"),
panel.background = element_rect(fill = "transparent"), # bg of the panel
plot.background = element_rect(fill = "transparent", color = NA), # bg of the plot
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text=element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank())
# export png file
overlay_file_wind <- "../data/fine_gridded_forecast_wind.png"
png(overlay_file_wind, width=image_size$width, height=image_size$height, units="px", bg="transparent")
p_wind
dev.off()
null device
1
img <- png::readPNG(overlay_file_wind)
plot.new()
rasterImage(img, 0,0,1,1)
在地形上叠加上述生成的地图, 气温分布以及风场流线图.
overlay_img <- png::readPNG(overlay_file)
overlay_img_temp <- png::readPNG(overlay_file_temp)
overlay_img_wind <- png::readPNG(overlay_file_wind)
# 2D plot with map overlay
shademat %>%
add_overlay(overlay_img, alphalayer = 0.95) %>%
add_overlay(overlay_img_temp, alphalayer = 0.95) %>%
add_overlay(overlay_img_wind, alphalayer = 0.95) %>%
add_shadow(raymat, max_darken = 0.4) %>%
add_shadow(ambmat, max_darken = 0.4) %>%
plot_map()
显示三维地形图.
zscale <- 10
rgl::clear3d()
shademat %>%
add_overlay(overlay_img, alphalayer = 0.95) %>%
add_overlay(overlay_img_temp, alphalayer = 0.6) %>%
add_overlay(overlay_img_wind, alphalayer = 0.8) %>%
add_shadow(raymat, max_darken = 0.4) %>%
add_shadow(ambmat, max_darken = 0.4) %>%
plot_3d(elev_matrix, zscale = zscale, windowsize = c(1000, 1000), background = "white",
water = FALSE, soliddepth = -150, wateralpha = 0,
theta = 25, phi = 30, zoom = 0.6, fov = 60)
label <- list(text="Olympic Stadium")
label$pos <- find_image_coordinates(long=116.39200, lat=39.99150, bbox=bbox,
image_width=dim(elev_matrix)[1], image_height=dim(elev_matrix)[2])
render_label(elev_matrix, x = label$pos$x, y = label$pos$y, z = 6000,
zscale = zscale, text = label$text, textsize=1.2, linewidth = 5, freetype=TRUE)
dist = distm(c(extent(elev_img)[1], extent(elev_img)[3]),
c(extent(elev_img)[1], extent(elev_img)[4]), fun = distHaversine)[,1]
render_scalebar(limits=round(dist/1000.,1),label_unit = 'km')
render_compass(position = "N",compass_radius=120, scale_distance = 1.2)
rgl::par3d(windowRect=c(-450, -100, 1000, 1100))
rgl::view3d(theta=20, phi=30, fov=45, zoom=0.6)
Sys.sleep(3)
render_snapshot('data/winter_olympic/output_image.png', title_text = "2022北京冬奥精细化预报",
title_size = 50, title_color = "black", title_font = "SimHei")
# rgl::rgl.close()